import numpy as np

def romberg_integration(a, b, f, epsilon=1e-5):
    from numerical_integration import composite_trapezoidal_fixed

    R = np.zeros([100, 100])

    # 初始化列表，包含 R00, R01 和 R11
    R[0][0] = composite_trapezoidal_fixed(a, b, f, 1)  
    R[0][1] = composite_trapezoidal_fixed(a, b, f, 2)  
    R[1][1] = (4 * R[0][1] - R[0][0]) / 3
    
    k = 1  
    while abs(R[k][k] - R[k - 1][k - 1]) > epsilon:
        k += 1
        
        R[k][0] = composite_trapezoidal_fixed(a, b, f, 2 ** k)
        
        # 外推过程 - 更新 R[k] 的内容
        for m in range(1, k + 1): # 下三角矩阵
            R[k][m] = (4**m * R[k][m - 1] - R[k - 1][m - 1]) / (4**m - 1)
    
    return R[k][k], 2 ** k
